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We present efficient quantum algorithms for simulating time-dependent Hamiltonian evolution of 
general input states using an oracular model of a quantum computer. Our algorithms use either 
constant or adaptively chosen time steps and are significant because they are the first to have 
time-complexities that are comparable to the best known methods for simulating time-independent 
Hamiltonian evolution, given appropriate smoothness criteria on the Hamiltonian are satisfied. We 
provide a thorough cost analysis of these algorithms that considers discretizion errors in both the 
time and the representation of the Hamiltonian. In addition, we provide the first upper bounds for 
the error in Lie- Trotter-Suzuki approximations to unitary evolution operators, that use adaptively 
chosen time steps. 



I. INTRODUCTION 

The original motivation for quantum computers stemmed from Feynman's famous conjecture that quantum com- 
puters could efficiently simulate quantum physical systems [l|, whereas there is no known way to achieve this with 
classical computers. This conjecture has spurred the construction of a number of quantum algorithms to efficiently 
simulate quantum systems under a Hamiltonian Q-Q- However, these algorithms are primarily for time-independent 
Hamiltonians. A simple extension to time-dependent Hamiltonians yields complexity scaling quadratically with the 
simulation time [l(| , a significant performance reduction over the near- linear scaling for the time-independent case [§1 . 
These issues can be resolved by generalizing the Lie-Trotter-Suzuki product formulae to apply in the time-dependent 
case. Such formulae have already been developed [Til [l2j . but have not yet been applied to quantum simulation algo- 
rithms. Here we explicitly show how these formulae can be used in quantum algorithms to simulate time-dependent 
Hamiltonians with complexity near-linear in the simulation time. We provide a number of improvements to further 
improve the efficiency, and a carefully accounting of the computational resources used in the simulation. 

Lloyd was the first to propose an explicit quantum algorithm for simulating Hamiltonian evolution |2J. This 
algorithm is for systems that are composed of subsystems of limited dimension, with a time-independent Hamiltonian 
consisting of a sum of interaction terms. The algorithm uses the Trotter formula to express the time evolution operator 
as a sequence of exponentials of these interaction Hamiltonians, which may be simulated efficiently. As a result, the 
complexity of the algorithm scales as 0(\\H || At) 2 , where At is the evolution time, and \\H\\ is spectral norm of the 
Hamiltonian. 

Aharonov and Ta-Shma @ and Childs Q extended these ideas to apply to Hamiltonians that are sparse, but 
have no evident tensor product structure. They use graph coloring techniques to decompose the Hamiltonian into 
a sum of one-sparse Hamiltonians, and use the Trotter formula and the Strang splitting respectively, to write 
the evolution operator as a sequence of one-sparse exponentials. The resulting sequence of exponentials can then 
be performed by a quantum computer. The use of higher-order splitting formula reduces the complexity of Child's 
algorithm to 0(\\H\\At) 3 ^ 2 , and it was conjectured that even higher-order formulae may lead to near-linear scaling pj. 

This h ypo thesis was verified by Berry, Ahokas, Cleve and Sanders (BACS) Q. They used Lie-Trotter-Suzuki 
formulae [14j j to generate arbitrarily high-order product formula approximations to the time-evolution operator, and 
gave an improved method for decomposing the Hamiltonian. The use of the Lie-Trotter-Suzuki formulae reduced the 
cost of their algorithm, causing it to scale as (||-ff|| At) 1+ °( 1 \ An alternative approach using quantum walks can yield 
scaling strictly linear in \\H\\At [H[2l[. 

High-order Trotter-like approximations for ordered operator exponentials are needed to extend the results of BACS 
to apply to the simulation of time-dependent Hamiltonian evolution. Such integrators were originally proposed by 
Suzuki pH . using a time-displacement superoperator. This method is made rigorous in Ref. [12], where sufficiency 
criteria for the applicability of the formulae, as well as bounds for the error, are provided. 

Here we explicitly apply these integrators to provide an algorithm for simulation of sparse time-dependent Hamil- 
tonians, and find that its complexity depends on the norms of H (t) and its derivatives. We show how adaptive time 
steps may be employed such that the complexity depends on average values of these norms, rather than the maximum 
values. This approach provides improved efficiency in situations where the norms have a sharp peak, or a finite 
number of discontinuities. For situations with singularities, we show how efficiency may be improved by adapting the 
order of the integrators. 
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We also improve the performance by specifying that the oracles that encode the Hamiltonian encode their outputs in 
polar form. Given this encoding, the one-sparse exponentials may be implemented via a simple circuit. In addition, we 
improve simulation efficiency by expressing the Hamiltonian as a sum in different bases. We quantify the performance 
of our scheme by considering the errors that occur in every step of the algorithm, including errors that occur because 
of discretization of the times used by our quantum oracles. We provide a unified presentation taking account of all 
these factors, as they interact in nontrivial ways. 



II. OUR APPROACH 



In this section, we provide a less technical explanation of the results and how they are obtained. Then we give the 
rigorous proofs in the following sections. The objective in a quantum simulation algorithm is to simulate evolution 
under the Schrodinger equation 



d_ 

df 



(i) 



where H{t) is the time-dependent Hamiltonian. That is, the initial state \tp(to)) is encoded in the qubits of the 
quantum computer, and we wish to obtain a state in the quantum computer encoding an approximation of the final 
state \ip{to + At)). A quantum computer simulation algorithm achieves this by applying an (encoded) approximation 
of the time-ordered exponential 



-i J H(u)du \ , 



(2) 



to the initial state in the quantum computer such that \ip(t + At)) = U(t + At, to) \ip(to)}. Given any e > 0, our 
goal is to obtain an approximation of the final state that is within trace distance e of the true state. This can be 
achieved Q if the approximation, U(t + At, to), satisfies 



C/(to + At,to)-C/(to + At,t ) <e, 



(3) 



with II • II defined to be the two- norm. 



A. Constant-Sized Time Step Simulation 



Our primary objective in this paper is to demonstrate a quantum simulation algorithm for time-dependent Hamilto- 
nian evolution that has a time-complexity that scales as At 1+ °W. The simplest approach to do so involves combining 
the sparse Hamiltonian decomposition scheme of BACS H, together with the higher-order integrators from Refs. 
[ill [T2I . BACS show that a Hamiltonian with sparseness parameter (the maximum number of nonzero elements in 
any nonzero row or column) of d may be decomposed into 6d 2 one-sparse Hamiltonians. Given that the state is 
encoded on n qubits, there is an additional factor of log* n to the number of queries required for the simulation. Here 
log* represents the iterated logarithm function, and increases extremely slowly with n. This factor arises because the 
decomposition requires 0(log* n) queries to the oracle for the Hamiltonian to perform the decomposition. 

The BACS decomposition technique may be used in the time-dependent case. However, one complication is that 
the sparseness can, in the completely general case, depend on time. That is, the nonzero elements at one time can be 
different to those at a different time. This would mean that the decomposition depends on the time, which makes the 
use of Lie-Trotter-Suzuki formulae problematic. To avoid that problem, we consider every matrix element that ever 
attains a nonzero value to be non-zero and use the BACS decomposition algorithm on those matrix elements. This 
allows us to directly apply their decomposition result. 

Reference [12| gives a result for exponentials of a general operator A(t). The result for Hamiltonian evolution 
follows by taking A(t) — iH(t). Then, for H(t) — Y^jLi Hj{t), by using a Lie- Trotter-Suzuki product formula that is 
accurate to order 2k fTll j . simulation error within e may be achieved using a number of one-sparse exponentials that 
scales as 

O (mk (y \ k (AAt) 1+1 ^/ e i/2M . (4) 
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Here A is an upper bound on the derivatives of the Hamiltonians such that 

/ \ V(p+i) 

A> (t)||j , (5) 

for t € [to, t + At] and p 6 {0, . . . , 2k} (25[. The notation with superscript (p) denotes repeated derivatives. An upper 
bound is used, rather than the exact maximum value, because the oracles that are used only give matrix elements of 
the Hamiltonian, not the norm. This is significant because methods for computing the norm of a matrix are often 
inefficient. However, it is often possible to place an upper bound on the norm, even if it is not possible to determine 
A exactly. 

We convert this result into a number of oracle queries by multiplying the number of exponentials by the number 
of oracle queries that are needed to simulate a one-sparse operator exponential. By doing so, we find that if the 
Hamiltonian is sufficiently smooth then the query complexity of the algorithm scales as (AAi) 1+ °W. 

More generally, we also consider the case where H (t) has discontinuous derivatives at a finite number of times. Such 
discontinuities are problematic because if a Lie-Trotter-Suzuki formula is used to integrate across such a discontinuity, 
then error bounds proved in [l2[ may not apply. In some cases this can be rectified by reducing the order of the 
integrator, but this strategy is not applicable if the Hamiltonian is not at least twice differentiable. Instead, we choose 
the time intervals to omit these points of discontinuity. In order to use this approach, it is necessary that the norm 
of the Hamiltonian is adequately bounded, because otherwise there could be significant evolution of the system very 
close to the point of discontinuity. Given this restriction it is possible to perform the simulation with complexity that 
is essentially unchanged. The full result, with the required conditions, is given in Corollary [6] 

It is important to note that the performance of our constant step size algorithm scales with the largest possible value 
of the norms of Hj and their derivatives. For some Hamiltonians, these values may only be large for a small fraction 
of the simulated evolution and so the algorithm may be inefficient. Using adaptive time steps, we can overcome this 
problem and demonstrate complexities that scale with the average values of the norms of Hj and their derivatives, 
given additional restrictions on the Hamiltonian are met. We discuss this approach below. 



B. Adaptive time steps 

The above scaling is the direct application of the results of Refs. @ and [l2j ■ We improve this scaling for problems 
where H(t) is badly behaved. For example, the matrix H(t) may be rapidly changing at some times and may be slowly 
varying at others. In such cases, using constant step size methods may be inefficient because overly conservative time 
steps will be taken during time intervals in which the Hamiltonian is comparatively easy to simulate. We address this 
by introducing adaptive time steps. When using adaptive time steps, instead of choosing each time step to be At/r, 
a sequence of times {t p } are chosen such that t < t\ < ■ ■ ■ < t r = t n + At. The size of the time intervals can be 
varied, as can the order of the product formula within each interval. 

We choose the duration of the time steps using a time-dependent function, T(i), that provides similar information 
to A, but at a specific time. We express this function as 

T (*)^ (eH p) w||J ' ( 6 ) 

for p G {0, . . . , 2k}. Throughout this work we use the notation for the average value, T(t&, t a ), defined by 



T(t b ,t a ):=- / T(t)dt. (7) 

tb — t a Jt a 

The goal is to choose the time steps adaptively such that the number of queries depends on the average value of Y(t) 
over the interval, rather than its maximum value (previously denoted A). The full result is given in Theorem [8] The 
basis of the method is to choose r time intervals to limit the error within each time interval to be no greater than e/r. 

To choose these time intervals appropriately, we need to know what the maximum value of T(t) is over a given 
interval, because the maximum value of T(t) dictates the simulation error over a short time step. Ideally one would 
want a method of choosing the duration of these steps that depends only on the value of T(t) at the beginning of the 
interval, in order to avoid needing to know the value of T(t) over the entire interval. We achieve this by requiring 
that the derivative of T(t) is appropriately bounded. A bound on the derivative is also necessary to obtain a result 
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depending on the average of T(t). This is because we can demonstrate that if the value of T(i) is approximately 
constant during each time step, then the complexity of the algorithm scales with the average value of T(t) over all 
time steps. We then require that the derivative of T(t) is bounded in order to guarantee this approximate constancy 
in the limit of short time steps. 

Even given the restriction on the derivative, it is unclear how to effectively choose the time intervals. The problem 
is that the time intervals are chosen to ensure that the error in each interval is no greater than e/r, but then the 
number of intervals will depend on how the intervals were chosen. To break this circular logic, we need a way of 
choosing an r g , such that when we choose the time intervals to have error no greater than e/r s , the total number of 
intervals is no larger than r g . The full technique is given in the proof of Theorem [8l The value of r g is chosen as in 
Eq. (|50p . The duration of each time step can then be calculated using just information about T(i) at the beginning 
of the time step, via an increment inversely proportional to T(t p ), as given in Eq. (j5l]). 



C. Resource Analysis 



Our cost analysis of these simulation algorithms focuses on the number of queries that are made to a pair of quantum 
oracles that provide information about the locations and values of the nonzero matrix elements of the Hamiltonian. 
We provide improvements to these oracles that enable us to improve the efficiency of the simulation, as well as to 
more precisely quantify the resource usage. First, we require that the output of the oracle for the values of the matrix 
is encoded using a qubit string as a complex number in polar form. The advantage of this is that the one-sparse 
Hamiltonian evolution can be applied using a simple circuit with qubit rotations proportional to the magnitude and 
phase of the matrix element. This is shown in Sec. IVIII1 and the explicit circuit is given in Fig. [3J Furthermore, the 
qubit rotations may be performed independently for each qubit of precision yielded by the oracle. 

To more precisely quantify the resource usage, we examine the number of qubits that the oracles need to provide, as 
well as the number of bits needed to represent the time. The oracles that are traditionally used in quantum simulation 
algorithms yield many-qubit approximations to the matrix elements in a single query. The fact that the qubits yielded 
by the oracle may be used independently motivates using oracles that output only one qubit per query. Doing so 
further reduces the number of qubits needed to simulate a one-sparse Hamiltonian evolution, because qubits that are 
not currently being used need not be stored. 

We find that the total number of qubits accessed for the positions of the nonzero matrix elements scales as n log* n 
(see Lemma [9]) , due to the need to access each of the n qubits for the position. This yields a factor of n increase in 
the apparent complexity over that if all qubits can be obtained in a single query. The number of qubits accessed for 
the values of the matrix elements is independent of n, but it does depend on other simulation parameters such as 
the error tolerance, the evolution time, the sparseness of the Hamiltonian, the norm of the derivative of H and k. It 
depends on all of these quantities logarithmically, except for k. The precise result for the number of oracle queries 
used is given in Lemma [1] 

We also consider a new source of error that is unique to the simulation of time-dependent Hamiltonians: the 
discretization of the time. Because the Hamiltonian varies with time, inaccuracy in the time will result in inaccuracy 
in the estimate of the Hamiltonian. This is potentially problematic for Lie- Trotter-Suzuki product formulae, which 
rely upon precise times in order to obtain higher-order scaling for the error. The higher-order scaling is not strictly 
obtained when the time is discretized, so it is necessary to show that the product formulae are not overly sensitive 
to error in the time so that we can use our discrete oracle in place of the continuous Hamiltonian. Another source 
of difficulty is in simulating systems with discontinuities. The problem is that the technique to avoid discontinuities 
requires choosing times arbitrarily close to, but on one side of, the discontinuity. With discretization of the time, the 
rounding may yield times on the wrong side of the discontinuity. 

In contrast to the output from the oracle, there is no need to use a coherent superposition of times, and the time 
may be regarded as a purely classical quantity at all stages in the calculation. This means that it is less challenging 
to provide the time to high accuracy than it is to obtain output from the oracle to high accuracy. Nonetheless, it 
is important for the reliability of the simulation that it is not unstable with inaccuracy in the time. We find that 
the simulation is stable with the time precision. The precision required depends on k, d, At and the maximum norm 
of the derivative of the Hamiltonian. It is logarithmic in all these quantities except for k; see Lemma [1] for the full 
result. The time precision also affects the result for simulating evolution with discontinuities in Corollary [5J That 
result holds provided the time discretization is no greater that the time between discontinuities (condition 4). 
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III. BACKGROUND AND DEFINITIONS 



Next we describe the technical background and definitions needed to understand our full results, which are given 
in the next section. For generality, we consider a Hamiltonian that is not sparse in any known basis, but is the sum 
of Hamiltonians that are each sparse in their own canonical basis. That is, 



M 



(8) 



where the set of operators {H' : 



i — ^ 



iNxN . 



[i= 1, . . . , M} is sparse when represented in their canonical bases. We 



express H'(t) = T^H^^T^, where H^(t) is sparse in the computational basis, and T M is a basis transformation that 
maps the computational basis to the canonical basis of H^t). 

To avoid complications due to the nonzero elements changing as a function of time, we consider only those elements 
which are nonzero at any time. We then take d to be an upper bound on the number of elements in any row that are 
nonzero at any time in the interval of interest: 



Definition 1. The set of operators {H^}, where H^ : R H- C 



NxN 



is d-sparse on S C R if for each fi there are at 



most d matrix elements in each row of H^(t) that attain a nonzero value for any t£ 5. 

Given a sparse Hamiltonian, BACS |8[ provide a method to decompose the Hamiltonian into a sum of 6d 2 one-sparse 
Hamiltonians, and express the evolution as a product of evolutions under each of these one-sparse Hamiltonians. The 
definition of sparseness we use here is compatible with that of BACS, and therefore each Hamiltonian H^{t) may be 
expressed as a sum of one-sparse Hamiltonians. That is, 



M 6d 2 

H(t) = Y,Y, T t H ^ T ^ 

[1=1 j=\ 



(9) 



where each H^j^t) is one-sparse. 

This decomposition is enabled by a quantum oracle that answers queries about the locations and values of the 
nonzero matrix elements of the Hamiltonian [6-9] . For the BACS decomposition, the matrix elements of each H^j (t) 
in ^ can be calculated using 0(log* n) queries to a quantum oracle for H^it). This oracle yields the locations and 
values of the nonzero matrix elements in a specified row of H^ to arbitrary precision In this case the cost of 
computing the matrix elements to sufficient precision is concealed within the definition of the oracle. To explicitly 
take account of the number of qubits that the oracle must provide, we separate it into two oracles that yield the 
positions and values of the nonzero matrix elements, and yield one qubit per query. See Sec. IV Al for the explicit form 
of the oracles. 

To express the evolution under the Hamiltonian as a product of evolutions under the one-sparse Hamiltonians, 
BACS use the product formulae of Suzuki 1 14] ]. Suzuki first proposed arbitrary-order product formulae for both time- 



IKIU 

independent and time-dependent cases [IJ, [14|. These product formulae are for operator exponentials, and are not 



limited to Hamiltonians. Subsequent work by the present authors showed that Suzuki's approximation method may 
be less accurate than expected if the operator does not vary sufficiently smoothly with time [l2j . That work provided 
upper bounds for the approximation error, given that the operator does satisfy a smoothness requirement. 

In this work, we use the result from Ref. [l2[ upper bounding the approximation error, with the operators obtained 
via the decomposition method of BACS. We therefore adopt the terminology "P-smooth" and "A-P-smooth" from 
Ref. These are formally stated below. 



Definition 2. The set of operators {Hj 



J 



= 1, 



7th 



1 c 



H<f\t) 



is finite on T. 



Definition 3. The set of operators {Hj 
Kp> (() ^> 



H 



j = I , . . . , m} is A-P-smooth on X C 
, for all t € T and p € {0, 1, • • • , P}. 



if, for each Hj, the quantity 



if {Hj} is P-smooth and A > 



The P-smooth requirement is needed in order to achieve an approximation of a given order, and the A-P-smooth 
requirement is needed to bound the error. In this work we also consider adaptive time steps, and then it is the upper 
bound on the derivatives as a function of time that is important. We therefore introduce the following definition. 



Definition 4. The set of operators {Hj : j 
R, if{Hj} is P-smooth and T(t) > f^™ 



1, . . . , m} is T -P -pointwise- smooth on the interval T C 
(p) (t)\\) , for all t el andpe {0,1,- •• ,P}. 



for T : 



i — ^ 
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In addition we adopt the terminology that a set of Hamiltonians is A-oo-smooth if the set is A-2/c-smooth for every 
k > 0. Similarly, we say that a set of Hamiltonians is T-cx)-pointwisc-smooth if it is Y-2fc-pomtwise-smooth for every 
k > 0. 

Provided {Hj} is A-2fc-smooth, for the Hamiltonian H{t) = YULi -ffijft) the evolution U(to + At, to) may be 
approximated via the integrator Uk, which is given iteratively via |12j | 

m 1 

Ux(t + At,t ) := Y[exp{-iHj(t +At/2)At/2} J| exp {-iHj(t + At/2)At/2} , 

j=l j=m 

U e {t Q + At, to) := Ut-i(to + At, to + (1 - Si )At)U^i(t + (1 - s^)At,t + (1 - 2*<)At) 

x U t -i(to + (1 - 2s*)At,t + 2s £ At)[/ £ _ 1 (£ + 2s e At,t + 8 t Ab)Ut-i(to + stAt,t Q ), (10) 

with S£ = 1/(4— 4 1 /( 2f ^ 1 )). This formula is implied by Suzuki's work [ll| . but is stated explicitly in [l2jj. The 
approximation error is 0((AAt) 2fe+1 ), so this formula is appropriate for short time intervals. For longer At, the 
evolution time may be divided into r subintervals, resulting in the approximation 

U(t + At,t )~f[U k (t +£—,t + (£-l)—} . (11) 

e=1 \ r r J 

The value of r is then chosen large enough such that the overall error is no greater than some allowable error, e. This 
choice of r scales as 0((AAi) 1+1 / 2fe /e 1 / 2fc ). More precisely, the error in the product formula will be no greater than e 
if we take 

r = \2e- l ' 2k (2fc(5/3) fc - 1 AAt) 1+1/2fc j , (12) 

provided that 

e< (9/10)(5/3) fc AA£. (13) 

This result is equivalent to Lemma 5 of [l2|, after eliminating Qk by using the inequalities in Eq. (A. 3) of that paper. 
The overall complexity of the simulation is then proportional to the value of r. 

IV. RESULTS 

This section formally presents our main results. The major result is an upper bound for the query complexity used 
to simulate time-dependent Hamiltonian evolution, using adaptive time steps and oracles that cost one query per 
yielded qubit. In order to quantify the complexity, the primary goal is to bound the error. In simulation schemes for 
time-dependent Hamiltonians there are three sources of error: 

1. integrator error from using Uk, 

2. error due to using a finite-bit representation of the time, and 

3. error due to using a discrctizcd representation of H^. 

To guarantee that the total error in the simulation is e, we ensure that the latter two errors sum to at most half the 
total error tolerance. The time steps are then chosen such that the contribution to the error from the integrators 
at most adds up to the remaining half. The following lemma, proved in Appendix |Al yields upper bounds for the 
number of bits of precision for these quantities that are needed to ensure that the roundoff errors are at most e/2. 

Lemma 1. Let {H^ : R H C 2 x2 ; (J, = 1, . . . , M} be a set of time- dependent Hermitian operators that is 2-smooth 
and d-sparse on I, which we take to be the union of disjoint closed intervals {Ij}. Furthermore define At := 
sup t . t, e x(t/ — U). The round-off error in simulating the product of the time- evolution operators that are generated 
by H{t) = ^^T^H^fyTfi over each Ij, using the integrator Uk is e/2 if 

1. the number of bits of precision used to represent the time, nt, and the number of qubits of precision used to 
represent the matrix elements, riH, satisfy 

•(max teI ^||a t i? A1 (t)||)(32fcMd 2 )(5/3) fe - 1 At 2, 



n t > 
n H >2 



log 2 
log 2 



32fcMd 2 (5/3) fc ~ 1 A max At 



6, (14) 



2. and for each j , the length of the subinterval Xj obeys 

\Zj\>At/Z 



(15) 



where A max is an upper bound on max ie i iA1 

Given these conditions on the number of bits of precision, we can then bound the number of queries to simulate 
the Hamiltonian using adaptive time steps as in the following theorem. 

Theorem 2. If {H^ : M i— > C 2 x2 ; fx = l,...,Af} is a set of time- dependent Hermitian operators that is T- 
2k-pointwise- smooth, d-sparse on I = [io^o + At], and nu satisfy (|14p , and i/iere exists K € M suc/i t/iat 
< -ftT 2 [Y(t)] 2 V ( £ I, i/iera /or any e G (0,1], £/ie evolution generated by H(t) = ^ZnT^H^t^T^ can be 
simulated within error e, and with the number of queries (denoted QueryCost) to our Hamiltonian oracles, satisfying 

QueryCost € O ( [n log* n + log(fcAf (5/3) fc d 2 A max At/e)] iV oxp ) , (16) 

( , , , Ad 2 T{t + M,t )M] 1+1,2k \ , s 

7V cxp € O ( Afd 2 fc(25/3) - K - J J , (17) 

where log* is £/ie iterated logarithm function and A max is given in Lemma\]\ The number of calls to {T^}, namely Nt, 
satisfies N T e O (N cxp / (3d 2 )) . 

This result is stated in asymptotic (big-O) notation wherein we take 

T(*o + At, to) At, A max Ai, max \\d t H^t)\\At 2 , M, d, n, e~\ k (18) 

to be our asymptotic parameters. Specifically, we say for two functions of these parameters, / and g, that / G 0(g) 
if there exists a constant a > such that |/| < a\g\ if all of these parameters are sufficiently large. 

If H(t) is the sum of sufficiently smooth terms and bounded on t E [to,oo), then k can be chosen such that (TI7j|) 
and (|17p scale nearly linearly with T(to + At, to) At. This observation is important because linear scaling is known to 
be a lower bound, and therefore our time scaling is nearly optimal 0, [22[. This observation is stated formally in the 
following corollary. 

Corollary 3. //, in addition to the requirements of Theorem^ {Hfj,} is T ' -oo-pointwise-smooth and d-sparse on [to, oo) 
then, 

QueryCost £ [nlog* n + n H ]Md 4 t(to + At, t )At(d 2 t(to + At, t )At/e)° (1) . (19) 

Before proceeding to show how to perform simulations with adaptive time steps, we first give the explicit scheme 
without adaptive time steps, and show how to take account of discontinuities in the Hamiltonian. 

V. SIMULATING TIME-DEPENDENT HAMILTONIANS 

In this section, we show how to simulate time-dependent Hamiltonian evolution on a quantum computer. In 
particular, we show that sparse time-dependent Hamiltonian evolution can be simulated efficiently provided that the 
Hamiltonian is at least piecewise twice-differentiable. We also show that if H (t) is sufficiently smooth, then the query 
complexity of our simulation scheme is comparable to the cost of the BACS algorithm for simulating time- independent 
Hamiltonian evolution. First we give the explicit description of the oracles used, then we give the precise result for 
the complexity of the simulation in terms of these oracles. 

A. Oracle Calls 

The time-dependent Hamiltonian over an Af-dimensional Hilbert space J^C can be represented by a matrix with 
elements H xy , for x the row number and y the column number. We consider a quantum oracle that can be queried 
to provide information about the locations and values of the nonzero matrix elements. For additional generality, in 
this work we assume that the Hamiltonian is in the form 



M 

H(t) = Y J T} l H li (t)T, l , 



(20) 
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where {H^} is d-sparse. This takes account of cases where the overall Hamiltonian is not sparse, but it may still be 
simulated efficiently using efficient basis transformations T M [3, H, H3| • We therefore require oracles to give the 
locations and values of the nonzero elements of the matrix representation of each (t) . 

Our first oracle provides the column numbers of the nonzero matrix elements in a given row of any H^. The function 
yields a requested bit of a particular entry in a list of d column numbers. This list contains the column number of 
every matrix element in a specified row of that attains a nonzero value. This function is Col, where Col(p, i, x, fj,) 
yields the p bit of the i th potentially nonzero matrix element in row x of H^. 

Our second oracle provides a requested matrix element of H^(t). This function yields a requested bit of a binary 
encoding of a given matrix element evaluated at a specified time. We denote this function as MatrixVal, and define 
MatrixVal(p, x, y, p, q) to yield the p th bit of a binary encoding of the matrix element [H fJ/ (t q )] x . To take account 
of discretization of the time, the time is specified by an integer q, which gives a time from a finite mesh {t q }, with 

t q =t + (q- 1/2) At/2 n \ (21) 

where n t is a positive integer and [to, to + At] contains the simulation time interval. We choose the matrix elements 
to be encoded in polar form, {H^(t)) xy = p(t) exp(i0(i)), where p and 4> are real numbers. For convenience, we also 
assume that p(t) is encoded as 

A^(f + f+- + fS). m 

where pj refers to the j th bit of p, and A max is an upper bound for max te i max (I= i. m ||-H/i(i)||max- 

The quantum oracles are unitary operations that give the values of these functions. That is, for classical inputs p, 
i and p and the quantum input | x) , 

QCol(p, i, p) \x) |0) = \x) \Co\{p, i, x, p)) . (23) 

This differs from Ref. Q, where i was given as a quantum input. In that work no superposition over the i was needed, 
so it does not change the analysis to give it as a classical input. Similarly, for classical inputs p, p and q and the 
quantum inputs \x) and \y), 

QMatrixVal(p, p, q) \x) \y) |0) - \x) \y) | MatrixVal (p, x, y, p, q)) . (24) 

In the following section, we present asymptotic estimates of the query complexity for simulating time-dependent 
Hamiltonian evolution using a sequence of approximations Uk, which are implemented on a quantum computer 
equipped with oracles QMatrixVal, QCol and {T^}. We will quantify the number of calls to QMatrixVal and 
QCol together as QueryCost, and quantify the number of calls to {!),} separately as iVr- 



B. Constant Timestep Simulation Method 



The simulation problem is as follows. Given the oracles QMatrixVal, QCol and the set of oracles {T M }, we wish to 
simulate the evolution generated by the Hamiltonian given in ([2"0")l . for {H^} A-2fc-smooth and d-sparse. In addition, 
the user is provided with an upper bound for the norm of each H^{t) and a similar upper bound for the norms of its 
derivatives. Our task is to provide an upper bound for the number of oracle queries that are needed to simulate the 
evolution generated by H(t) within error e. These upper bounds are given in the following lemma. 

Lemma 4. If {H^ : 1h> C 2 x2 ; p = 1, . . . , M} is a set of time- dependent Hermitian operators that is A-2k-smooth 
and d-sparse on [to, to + At], then for any e £ (0,1] the evolution generated by H(t) = T^H^^T^ can be simulated 
with the error due to the integrator bounded by e/2, and a number of queries to QCol and QMatrixVal satisfying 



QueryCost < 12CMd z 5 



2rfe-l 



2Akd 2 AAt 



6d 2 AAt \ 1/2k ' 
~Tj2T) 



(25) 



where e := min{e, 18(5/3) l d 2 KAt\, C is the number of oracle calls needed to simulate a one-sparse Hamiltonian. 
The number of queries to {T^}, Nt, is bounded above by QueryCost/ (3Cd 2 ) . 

Proof. Our approach is to express H(t) as a sum of 6Md 2 one-sparse Hamiltonians and use (|11[) to approximate 
U{to + At, t ) with a sequence of these one-sparse operator exponentials. The number of oracle calls in the simulation 
is then obtained by multiplying the number of one-sparse exponentials in our approximation by the number of oracle 
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calls to simulate one-sparse Hamiltonian evolution. That number is just given as C here, and an upper bound will be 
placed on it in Lemma [§] 

Each d-sparse Hamiltonian may be decomposed into a sum of 6d 2 one-sparse Hamiltonians H^j by using the 
BACS decomposition scheme. As each matrix element of is uniquely assigned to a one-sparse Hamiltonian H t 



in the BACS decomposition [8j, we have that \\H^-{t)\\ < for any non-negative integer p < 2k. Using 

Definition O if {H^} is A-2fc-smooth, then the set of Hamiltonians {H^ j} is 6d 2 A-2fc-smooth. Using Eq. (fT2")l and 

if e/2 < (9/10)(5/3) fc 6(i 2 AAt, then the number of exponentials needed to simulate 



(p), 



!>•} 



the fact that N. 



cxp 



2m5 



Hamiltonian evolution, using constant-sized time steps and within error e/2, is bounded above by 



N cxp < 2m5 



fe-i 



2[2fc(5/3) fe - 1 6d 2 AAi] 
(e/2)!/2fc 



1 + 1/2A: 



(26) 



Note that we have replaced e with e/2, because we require error bounded by e/2 here. To ensure that condition CE 
holds, we then replace e by e, which ensures that this condition holds and that the error is no greater than e/2. 
Using to = 6Md 2 and simplifying gives 



A^oxp < 12Md 2 5 k - 1 



2Akd 2 AAt - 



6d 2 AAt\ 



l/2k' 



(27) 



The number of oracle queries that are needed in our simulation is QueryCost = CN exp , which gives Eq. (J25J). 

Finally, we verify the claim that the number of basis transformations is bounded above by N CKp / (3d 2 ) by counting the 
number of basis transformations that result from using the Lie- Trotter-Suzuki formula. As the BACS decomposition 
method expresses a d-sparsc Hamiltonian as a sum of 6d 2 one-sparse Hamiltonians, and {H^} is d-sparse, it follows 
that the Hamiltonian can be expressed as 



M I 6<f 



(28) 



where H^(t) = Y^j=i ^/*,j(*)> an d {H^j} is one-sparse. We can then use (fT0|) to show that Ui(t a + r, t ) becomes 



M I 6d' z 

II T l ( II exp(-iff M>i (*o + r/2)r/2) ] T p 



f[ TU JJ exp(-iF M -(to + r/2)r/2) ] T, 



(29) 



Equation (|29|) has only AM basis transformations, but \2Md 2 one-sparse operator exponentials. Because Uk is a 
product of 5 fe_1 such approximations, there are at most 3d 2 basis transformations per one-sparse operator exponential 
in Uk- Because the approximation to U in our decomposition is a product of Uk [12j |. 



N T < N cxp /(3d 2 ), 



which implies that Nt < QueryCost/ (3Cd 2 ) via this method. 



(30) 
□ 



This Lemma shows how the complexity scales when the Hamiltonian is permitted to be a sum of terms that are 
each individually sparse in different bases. This result only considers the error in integrator, and does not consider 
the contribution of the round-off error that occurs due to discrctizing both time and the matrix elements of . The 
following Theorem gives values of the precision that are sufficient to ensure this round-off error is e/2, implying that 
the total error is at most e. 



Theorem 5. // {H^ : E H> 



;/i = l,... , M } is a set of time- dependent Hermitian operators that is d-sparse 



and A-2k-smooth on [to, to + At], and n t and nn satisfy (|14p . then for any e € (0,1], the evolution generated by 
H(t) — ^2 T^H^fyTft can be simulated within error e, while using a number of oracle queries to QMatrixVal 
and QCol, QueryCost, that are bounded above by 



QueryCost < YlCMdTh 



2rfe-l 



24kd 2 AAt 



6d 2 AAt\ 1/2k ' 



e/2 J 



(31) 



where e := min {e,18(5/3) fe - 1 d 2 AAt} and the number of queries to {T^}, Nt, obeys Nt < QueryCost/ (3Cd). 
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Proof. The error in our simulation scheme arises from two sources: the discretization error, and the error due to the 
integrator. By requiring that n t and tih satisfy (|14[) , Lemma [T] implies that the error due to the discretization is no 
more than e/2. Note that, because I is just a single time interval, the At here corresponds to the At in Lemma Q] and 
the condition (|15l) is automatically satisfied. Using Lemma HI we find that the simulation can be performed such that 
the error in the integrator is no greater than e/2, and the query complexity of the simulation satisfies the inequalities 
in Eqs. (|25|) and (|30[) . Because e < e, using the triangle inequality shows that the total error is no greater than e. □ 



In some cases where {H^} is not smooth, the Hamiltonian evolution can be simulated more efficiently by deleting a 
neighborhood from [to, to + At] around each point where the derivatives of {H^} diverge, and using the integrator Uk 
to approximate the time-evolution in the remainder of the interval. The query complexity for simulating Hamiltonian 
evolution by this method is given by the following corollary. 



- 2 ; /i = 1, . . . , M} be a set of time- dependent Hermitian operators that is A-2k- 

< tx, < + At, with the additional conditions 



Corollary 6. Let {H^ : R i-> 

smooth on I = (to, to + At) \ {t%, . . . , t^}, where to < t\ < 
1- 3 -ff max € K : H max > max i6 [ t0! j 0+At ] ||i/(t)||, 

2. < e < min{l,27(5/3) A; - 1 d 2 AAt} 7 

3. n t and nu satisfy (|14p . and 

4. At/2 nt < min^ =0 ,...,z,(^+i - t t ) with t L+1 := t + At. 

Then the query complexity for simulating evolution generated by H(t) = 'Yli^T^H ^(t)T^ within an error of e is 

k / 6d 2 AAt\ 1//2fe 



QueryCost < UCMd'h 



2rk-l 



(L + l) + 24fcd 2 AAt 



(e/3) ; 



(32) 



where C is the number of oracle calls needed to simulate a one-sparse Hamiltonian, and < QueryCost/ (3Cd 2 ) . 

Proof. We remove ^-neighborhoods around each tf and simulate evolution on the remaining time interval. We then 
exploit the fact that H(t) is bounded to estimate the error incurred by omitting the evolution around those points. 
We choose a value of S satisfying 



< S < min 



min (t l+x - t t ) - At/2 nt 

l=Q,...,L 



1 



' 2H„ 



log 



1 



e/6 



where i^+i := + At. The evolution operator is 



U(t + At, t Q ) = U(t L+1 ,t L+1 - 5) 



[I U(t e+1 -5,U + 5)U(t e + 5,ti-S))x Ufa -5,t + S)U(t + 6, t ) 



t=L 



l[u(t e+1 -5,t e + S). 



(33) 

(34) 
(35) 



For e < 27(5/3) fe - 1 d 2 AAt 



2e(t t+1 -t t - 26)/ (3 At) < 18(5/3) k ~ 1 d 2 A(t i+1 -tt- 26). (36) 

• )J can be simulated with integrator error 



With this restriction, using Lemma 2] each of the L + 1 evolutions in 
bounded above by e(tt+i — tt — 26)/ (3 At) using no more than 



2efc-l 



24fcd 2 A(^+i -tt- 26) 



\3 



6d 2 AAt 



l/2fc' 



(37) 



oracle queries. Using Eq. (4.69) of Ref. [T5(, which states that for unitary operators, 



U u i-U v * 



3 



v \l 



(38) 
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the total error is bounded above by 

^e{t l+l -t l -28)/{3At)<e/3. (39) 

e 

Then, summing ()37[) over I gives inequality (|32j) as an upper bound for the number of oracle queries made. The 
additional factor of (L + 1) in Eq. (|32|) is to take account of the ceiling function in Eq. (|37| . The bound for the 
number of basis transformations is obtained by using Nt — N cxp /3d 2 , and by summing over all L + 1 subintervals. 

To bound the overall error, we need to bound the error due to approximating (f34|) with (|35|) . Because H(t) is 
bounded on [to, to + At], and H max > max te [ tn .t +At] \\H(t)\\, the unitary evolution over each tt for I = 1,...,L 
satisfies 

\\U(tt + 5,t t -S)-l\\ <e 2H ™ s -l. (40) 

For £ — and L + l we have 

||£/(t + <Mo)-l|| <e H — 4 -l, 
\\U(t L+1 ,t L+1 -8)- 1|| < e ff — 5 - 1. (41) 

These errors can be made suitably small by using the restriction on 8 specified in Eq. (1331) . The first expression in 
the braces in Q33p ensures that the inequality 8 < min^ = o,... ~ te)/2 is satisfied. The second ensures that the 

error in approximating the evolution about each of the tt by 1 is bounded above by e/[6(L + 2)]. Using Eq. (|38[) . the 
total error in approximating (IMl) with (l3"5j) is bounded above by e/6. Combining this with the bound on the error 
of e/3 for the integrators over the time intervals \tg + 8, tt + \ — 8], the total error due to omitting the 8- neighborhoods 
from the simulation and using the Lie- Trotter-Suzuki formula is bounded above by e/2. 

We now use Lemma Q] to ensure that the roundoff error is also bounded above by e/2. The definition of At used 
in that Lemma gives the At used in this corollary, and so may be used in the restrictions without change. The 
restriction At/2 nt < min^(i£ + i — tt) and the choice of 8 in ([53")) ensure that the restriction (fTaj) of Lemma [1] holds. We 
have also required that n t and satisfy (I14[) , so all conditions required for Lemma [T] hold, and the round-off error 
may be bounded by e/2. As the error in the integrator has also been bounded by e/2, the total error is no greater 
than e. □ 

We have shown in this section that time-dependent Hamiltonian evolutions can be simulated by using the product 
formula approach to simulation, even if there are discontinuities in the Hamiltonian. If the Hamiltonian is sufficiently 
smooth, then these simulations achieve the same near-linear scaling as the BACS simulation achieves. In the next 
section we improve upon these results by presenting a method that uses adaptive time steps. 



VI. ADAPTIVE DECOMPOSITION SCHEME 



We saw in the previous section that if we use time steps that have constant size in our simulation algorithm, then 
the number of oracle calls used to simulate Hamiltonian evolution depends on the largest values of the norms of 
H(t) and its derivatives. For Hamiltonians whose time-dependence is sharply peaked, the maximum values of these 
quantities can be quite large relative to their time-averages. It is natural to ask if the complexity can be made to 
depend on the average values, rather than the maximum values, by using adaptive time steps. In this section we show 
that this is indeed possible. 

As discussed in Sec. Ill B[ this result is nontrivial because of the interdependence of the different parameters. We 
choose a sequence of times {t p } such that to < t\ < ■ ■ ■ < t r = to + At. These times are selected such that the error 
from using Uk is bounded above by e/(2r) for each interval. Given a function T, such that {H^} is T-2/c-pointwise- 
smooth on [to, to + At], the size of the interval will be inversely proportional to the maximum value of T(t) in that 
interval. The exact result is given in the following Lemma. 

Lemma 7. Let {H^ : /i = 1, ...M}, where : R i— > C 2 x2 , be a set of time- dependent Hermitian operators 
that is T-2k-pointwise-smooth and d-sparse on [to, to + At], and {t\, . . . ,t r } be a set of r moments in time such 
that t < t\ < ■ ■ ■ < t r = to + At. If 

/ \ f /_\l/(2*+l) 

U^/^^-^- ^W- 1 ' (42) 
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where e G (0, 1], then 



U(t + At, i ) -Y[U k (t p ,t p -x] 

P =i 



< e/2. 



Proof. Using Theorem 3 and Eq. (A. 3) of Ref. Qjl, the error in each interval is bound above by 

\\U(t p + At p , t p ) - U k (t p + At p , t p )\\ < 2[2fc(5/3) fe - 1 AAt p ] 2fe+1 , 
where At p :— i p+1 — t p , provided that 

4fcV2(5/3) fc ~ 1 AAi p < 3/2, 



(43) 



(44) 



(45) 



for a set of operators that is A-2fc-smooth. As A is bounded above by 6c? 2 max tG [ tp4p+1 ] T(t), the upper bound in (|4"4")l 
becomes 



|| U(t p + At p , t p ) - U k (t p + At p , t p )\\ < 2 



-i 2k+l 



12d 2 fc(5/3) 



k-l 



max T(t)At p 

t£[tp,tp+i] 



(46) 



Because e < 1 and r > 1, the condition (|4"2")l implies that (1431) is satisfied, and therefore that Eq. (|4"6"]l holds. 
Using Eq. ((42j) in Eq. ([46]) gives 



Then, using Eq. 



||L7-(t,, + A*p,tp)-E/k(tp + Afp,tp)|| <e/(2r). (47) 
the error due to using r Lie- Trotter-Suzuki integrators is at most e/2, hence yielding (|4"3")l . □ 



There are two challenges in using this result as a method for choosing the time steps. First, Eq. (l4"2l depends on 
the maximum value of T(t) over the interval, which may be difficult to find, particularly because it depends on the 
size of the time interval. Ideally, it should depend only on the value at time t p , in order to give a simple method 
of determining the time interval. Second, Eq. (I42[) depends on the number of time steps r, which is not known in 
advance. It is only known once all time intervals have been determined, but we wish to determine these via Eq. (|42l) . 

To break these interdependencies, we first need a bound on the derivative of T(t), so we can bound the value of T(t) 
on the interval by the value at time t p . It is not obvious what bound should be taken, but it can be shown that if Tp 
is the smallest possible function such that {H^} is T-P-pointwise-smooth, then it will satisfy \dtTp(t)\ < [Yp + i(t)] 2 
(see Appendix [Bj . In many cases we can expect that there exists a constant K such that Tp + i(t) < KTp(t). This 
motivates taking the condition |<9 t T(t)| < K 2 [T(t)] 2 . 

Second, we provide a quantity r ff , given in Eq. (|50p . that is an upper bound on the number of time steps. If the 
time steps are chosen according to Eq. (|4*21 with r g instead of r, then Eq. (|4*2"j) must still hold. Using this approach, 
we obtain the following Theorem. 



Theorem 8. // {Hp 



i — ^ 



-.2"x2" 



; /i = 1, . . . , M} is a set of time- dependent Hermitian operators that is T- 



2k-pointwise- smooth, d-sparse on I = [to, to + At], nt, and nn satisfy (|14l) . and there exists K G K such that 
|<9tT(i)| < K 2 [T(t)} 2 V t € I, then for any e G (0,1] 7 the evolution generated by H(t) = ^2pT^H^(t)T^ can be 
simulated within error e, with 



QueryCost < YlCMd 2 h 



1rk-\ 



[24d 2 fc(5/3) fc - 1 T(t + At, t ) At] 
{e/4y/ 2k 



l + l/2fe 



3AT 2 T(t + At,t )At + l 



(48) 



where log* is the iterated logarithm function, A max is given in Lemma [U C is the number of oracle calls needed to 
simulate a one-sparse Hamiltonian, and < QueryCost/ (3Cd 2 ). 



Proof. For this proof, we choose a set of r times {ti, . . . , t r }, such that t < t\ < ■ ■ ■ < t r = t + At. We define 

[24d 2 fc(5/3) fc - 1 T(t + At,t )At] 1+1 /2fc 



.4 



(49) 



and 



\A + 3K 2 T(t + At, t ) At + 1] . 



(50) 
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We choose the times via the recurrence relation, for p = 0, . . . , r — 1, 



where 



Y := 



T(t p ) Y + K 2 ! 

24d 2 fc(5/3) fc - 1 
( e / r y/(2k+i) ■ 



(51) 



(52) 



We will show shortly that (f5"Tj) implies that P2"f is held, thereby guaranteeing that the integrator error is bounded 
above by e/2. In order to prove this we first need to find an upper bound for T(t + St), for small values of St > 0, in 
terms of T(f). 

We find this bound by solving the differential equations given by T saturating the inequality |StT(f)| < K 2 [T(t)] 2 . 
In particular, for St > 0, we have 



T(t) 



1 + J ftT 2 T(i)(5t 



< T(t + <5t) < 



T(t) 



1 - K 2 T(t)St' 



(53) 



By using the inequality on the left, and the fact that the minimum of a function is less than or equal to its average 
value, we find 



[tp+1 > p) ~ l + K 2 T(t p )At p 



This, together with recurrence relation (J5TJ), yields 

1 <T(t p+1 ,t p )At p (Y + 2K 2 ). 
Summing over the first r — 1 subintervals then gives 

r-2 
p=0 

= (t r -l - to)T(*r- 1> t )(y + 2if 2 ). 



(54) 



(55) 



(56) 



By solving the above equation for r, and using the fact that (i r _i — to)T(i r _i,fo) < AtT(to + At, to), we have that 



r < AtT(t + At, t ) (Y + 2K 2 ) + 1. 
As r is an integer, it is also bounded above by 

r < [AtT(t + At, t Q ) (Y + 2K 2 ) + l\ . 

This also implies that for any 7 > 0, 

r < \ AiT (to + At, to) (Y + 2K 2 ) + 7] . 
If we take 7 = 1/3 in (|59p and use the definitions of A and Y, we find 



(57) 



(58) 



(59) 



r < 



< 



< 



A 2fc/(2fe+l) r l/(2fc+l) + 2K 2 Y{tQ + AL to)At + 1/3 " 

A[l + (r g /A - l)] 1 /(2fc+i) + 2K 2 T(t + At, t )At + 1/3 



A 



2k 



-(r g -A)+ 2K z T(to + At, to) At + 1/3 



A+\ (3K 2 T(t + At, to) At + 2) + 2K 2 T(t + At, to) At + 1/3 



< \A + 3K 2 T(t + At, f Q )Af + 1] 



(60) 
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FIG. 1. This figure shows a comparison between 
the number of exponentials used in our adap- 
tive and constant step size methods for H(t) = 
exp(-(t- l)7a 2 )/(a v / ¥)l, t G [0,2], e = 1(T 4 
and k — 2. 



FIG. 2. This plot shows the number of exponen- 
tials found by choosing k adaptively for H(t) = 
t 5 sin(l/i) exp(— t)l, and e = 10~ 4 for one and two 
subintervals with k = 1 and k = 2 respectively. 



The inequality on the right of Eq. (f55[) gives 



max T(t) < T }l p ) . t . (61) 

te[t p ,t p+1 ] w ~ 1 - K 2 T{t p )At p y ' 



This, together with the recurrence relation (|5Tj) . gives 



max T(t)At p < l/Y. (62) 

Because r < r s , the condition (|42|) of Lemma [7] is satisfied, and the error from the integrator is no more than e/2. 
This implies that the total error is bounded above by e if nt and nn are chosen as in Lemma [TJ 

We can place an upper bound on the number of exponentials used in the approximation by multiplying r by the 
number of exponentials in each A: th -order Lie- Trotter-Suzuki approximation. We use r g as an upper bound for r, note 
that there are 12Md 2 5 k ~ 1 one-sparse operator exponentials in each of the fc th -order Lie- Trotter Suzuki approximations, 
and use Query Cost = CN exp , giving Eq. (|35]) . The bound on iVr follows from the fact that Query Cost = CN exp , and 
N T < N cxp /3d 2 . ' □ 

The result of Theorem |8] shows that adaptive time steps can be used to dramatically reduce the complexity of simu- 
lating certain time-dependent Hamiltonian evolutions. This improvement stems from the fact that the performance of 
the constant time step algorithm depends on the largest possible value of T, whereas the performance of the adaptive 
version scales with the average value of T. In some cases, such as an example in the following section, the average 
can be much smaller than the largest value of T, leading to a substantial difference in performance. 



VII. EXAMPLES 



We now examine the performance of our adaptive time step decomposition method for a pair of examples. In 
the first example, we examine a Hamiltonian that is a Gaussian approximation to a Dirac-delta function. In our 
second example, we examine a non-analytic Hamiltonian, and show that we can substantially reduce the number of 
exponentials used in some simulations by choosing k adaptively as well. 

The Hamiltonian that we choose in our first example is, for a > 0, 

H(t) = exp(-(f - if /a 2 )/{a^)t. (63) 

We compare the cost by plotting -/V oxp as a function of a for both methods. We choose our approximations to be 
2 nd -order Lie- Trotter-Suzuki formulae (k — 2) with e = 10 -3 . For our adaptive method, rather than choosing the 
time steps using the upper bound on r g as in Theorem [51 we use an iterative process to find the number of steps. 
That is, we choose each step using (|4"2"j) with an initial guess of r — r g , then find a sequence of steps via (|4"2")l . and 
then count the number of steps to find a new guess for r. Iterating this process gives the solution for r. 
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Fig. [T] shows that the adaptive method results in a value of iV cxp that is approximately constant in a. In comparison, 
the number of exponentials used in the constant step size case diverges as a approaches zero. This divergence occurs 
because lim a ^o m ax tg m2] 1*4 (i) = oo, whereas the adaptive method yields a nearly constant value of N exp , because 

Jp 2 Ti{t)di only weakly depends on a. This shows that it can be advantageous to choose the step size adaptively if 
H (u) varies substantially over the interval of simulation. 
The Hamiltonian that we use in our second example is 

H(t) = t 5 sin(l/t) exp(-t)l. (64) 

The second-order Lie- Trotter-Suzuki formula should not be used as an approximation to J7(At, 0) [l^ |. because the 
third derivative of H(t) diverges near t = 0. However, XJi can be used to approximate the time-evolution on any closed 
interval that excludes t = 0. A natural way to handle this is to divide the approximation into two subintervals, and 
use different values of k to approximate the evolution within each subinterval. The evolution in the first subinterval, 
u G [0,f'], is approximated using U\, whereas the evolution in the remaining subinterval is approximated using L^- 
We reduce the number of exponentials that are used in our simulation schemes for any fixed value of t' by applying 
our adaptive step size algorithm in each subinterval. We then vary t' using a gradient search method to change the 
size of each subinterval to further reduce the number of exponentials that are used in the simulation. 

We show in Fig.[5]that choosing k adaptively can lead to reductions in the number of exponentials that are needed to 
approximate the time-evolution operator. Therefore, when simulating the evolution generated by Hamiltonians with 
singularities, it may be more efficient to use lower-order formulae near the singularity, and higher-order formulas further 
away from it. This method may have applications in situations where high-order integrators fail to provide the scaling 
expected for singular differential equations, such as those that occur in Coulomb problems [23|. 

VIII. SIMULATING ONE-SPARSE HAMILTONIANS 

In Sections [V] and I VII we specified the cost of simulating the evolution as a function of C, which is the number of 
oracle calls that are needed to simulate a one-sparse Hamiltonian. In this section, we provide an upper bound for C 
and discuss how this cost relates to those obtained using other oracle definitions. We then use this bound for C in 
concert with our results from Sec. IVII to prove Theorem [5] and Corollary [3J 

Lemma 9. Let {H^ : fj, = 1, . , . , AI}, where : R n- C 2 x2 , be a set of time- dependent Hermitian operators that 
is 2k-smooth on [to, to + Af], and let {H^j} represent the set of one-sparse Hamiltonians that result from applying 
the BACS decomposition algorithm to each {H^}. The query complexity, C, of simulating exp {— iH^ t3 (r)St} for any 
j and [r, r + St] C [to, to + At] using the oracles QCol and QMatrixVal and nn bits of precision to represent the 
matrix elements, is bounded above by 

C <4n(z n + 2) + 3n H , (65) 

where z n is the number of times the mapping z i— > |~21og 2 (z)] must be iterated, starting at n, before reaching a value 
that is at most 6. 

To prove this Lemma, we take advantage of the polar decomposition of the matrix elements of the Hamiltonians. 
This enables a remarkably efficient simulation of the evolution under one-sparse operator exponentials using a sequence 
of rotations that are controlled by only one qubit at a time. In particular, the result is as in the following Lemma. 

Lemma 10. Given that the assumptions of Lemma\Qare met, and if nn bits of precision are used to represent each 
of the matrix elements of Hj, then exp(— iH^ j (r)5t) can be simulated using one qubit to store these matrix elements 
and 3nn queries to the oracle QMatrixVal. 

Proof of Lemma \l(A The decomposition scheme of BACS Q takes a given row number x, and determines if there will 
be a matrix element in row x assigned to H^ j. If there is, then the boolean function £(x) is set to 1; otherwise it is 
zero. If there is a matrix element, then the column number is determined, and M x and m x are set equal to the row 
and column numbers, such that M x > m x . The decomposition method also generates a three-bit string v{x), which 
we will not otherwise use in the simulation scheme. The BACS decomposition method will therefore transform the 
initial state 1-0) according to 

|0) = J2 «x \x, 82 "+ 4 ) ]T a x \x, m x , M x , v(x),£(x)) . (66) 

X X 
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FIG. 3. This circuit simulates exp(—iH tJ , : j(tp)8tp) for the one-sparse Hamiltonian H^j, given an input state of the form of 
the LHS of (167)1 . Here the variable 4> = Arg([// M ,_, (t p )] mx ,M x ) and a = 2| [H^j (t p )]m x ,M x | Ai p . Here we also use rectangles to 
represent R z rotations by a fixed angle, triangles represent R x rotations and the pentagon represents a R y rotation. These 
rotations can be enacted by querying the oracle QMatrixVal and performing controlled rotations on the output. 

Then given this transformed state, our next goal is to perform a mapping between the subspace span{|ma;) , |-Mx)} 
and an ancilla qubit space. The purpose of this is to map this two-dimensional subspace onto one that can be evolved 
using single-qubit operations. This transformation is 

( x = M x = m x , \m x © M Xl M Xl v(m x ),£(x), 1, 1) 
\x,mx,M x ,i>{m x ),£{m x ),0,Q) H> ^a x |0 n ) ® I x = M x ^m x , \m x © M x , M x , v(m x ), £{x), 1, 0) . (67) 
x x [x = m x ^M x , \m x © M x> M x ,v(m x ),£(x),Q,0) 

The second last qubit is the one that encodes the two-dimensional subspace; it is |1) if x = M x , and |0) if x = m x . 
The last qubit is used to indicate if a; is a member of a one-dimensional subspace. 

Given a state in the form of the RHS of (|67p . we can simulate the evolution of |a:) by evolving the second last 
ancilla qubit, whose state is logically equivalent to \x). To do so, we must know whether x is in a one- or two- 
dimensional irreducible subspace. Specifically, the evolution operator takes one of two possible forms on the subspace 
span(|m x ) , \M X )). It performs the transformation 

ajf, ^ «m x exp(-i[H" M (i p )] majMa; At) (68) 



«M 7 



exp 





[-ff M (ip)]m x M x 







Atj 









(69) 



if the subspace is one- or two-dimensional, respectively. 

We can write the two-dimensional rotation as a sequence of Pauli-z and -y rotations using standard decomposition 
techniques 15]. Given i? TOx! j\/ x = pexp(i^>), the resulting decomposition of the two-dimensional transformation is 



exp 







m x M x 







At r 



R z (-7r/2)R z (- ( j ) )R y (2pAt p )RMR z ^/ l 2)- 



(70) 



The sequence of exponentials in (|70l) can be implemented using a sequence of Pauli rotations that are controlled by 
the qubits that encode the matrix element £f mxi M x - 
The one-dimensional subspace is evolved according to 



(71) 



exp -i[if M (i p )] xa; A*p | a;) . 



Because we have encoded \x) in the one-dimensional case to be |1) in (|67[) . the one-dimensional transformation can 
be expressed as R z ( — 2[H ll (t p )\M x M lc At p ). This rotation can also be written as, 



R z (- 2[H tl {t p )] Mx M x At p ) = R z (-<p)R x (-7r/2)R y (2pAt p )R x (7r/2)R z {4 > ), 



(72) 
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which allows us to write the one-dimensional rotation in a form that is similar to the two-dimensional rotation. 

We present a circuit in Fig. [3J that enacts both the transformation in (|67p and also (f72"j) and (J7UJ) coherently on 
each subspace. We combine the rotations for both the one- and two-dimensional cases together in a single sequence 
of rotations, by making the tt/2 rotations controlled by the last qubit. In addition, we allow the rotations to be 
controlled by £(ir), so no rotations are performed if no matrix element has been assigned to row x of H^j. 

The rotations R y (a), where a = 2pAt p , and R z {4>) in Fig. [3J can each be implemented by calling the oracle 
QMatrixVal nujl times, provided uh is even, and equal numbers of bits are used to encode the modulus and phase 
of the matrix element. Since there are three rotations of this form, 3n#/2 oracle calls can be used to enact them. 
However, in order to re- use the ancilla bits that record these values, we need to call the oracle another 3nn/2 times. 
Therefore, the total number of calls made to QMatrixVal is bounded above by 3n# . □ 

Note that, in Fig. [3] we can access each qubit of the oracles independently, without needing to store the other qubits. 
This is because, for each controlled operation (such as R z ((f>)), we can call the oracle for one qubit of precision, perform 
the rotation for that qubit, then call the oracle again to erase the value, before calling the oracle for the next qubit. 
Now that we have proven this Lemma, the proof of Lemma [S] is simple. 

Proof of Lemma\Q The one-sparse matrix exponentials may be performed by using the BACS decomposition tech- 
nique, then performing the rotations as described in the proof of Lemma [TU1 then inverting the BACS decomposition 
technique to restore the ancilla qubits to their original states. The BACS decomposition technique uses 2(z„ + 2) 
queries_to their oracle to identify the irreducible subspace that a basis state is in, and store this information in a qubit 
string Using an oracle that only provides one qubit at a time, the number of oracle calls is multiplied by a factor 
of n. Another factor of 2 is obtained because the BACS decomposition technique is inverted, yielding a total number 
of queries of An(z n + 2). Using Lemma ITUI the rotations can be performed using 3nn calls to QMatrixVal, so the 
total number of oracle calls needed is bounded as in Eq. (|65p . □ 

Using Lemma it is now straightforward to prove Theorem [5J which is our main result in the paper. 

Proof of Theorem, [H The proof follows directly by substituting the result of Lemma |H] into those of Theorem [FJ while 
noting that the second and third term in (|48l) are asymptotically subdominant. □ 

As discussed in Ref. [l2j for the case of constant time steps, if the Hamiltonian is sufficiently smooth then we can 
choose k to increase with At, so that the complexity scales close to linearly in A At. Here we obtain a similar result 
for the case where the time steps are chosen adaptive ly. Theorem [2] provides a guide to choose an optimal value of k, 
which then enables us to prove Corollary [3J 

Proof of Corollary^ As {H^} is A-oo-smooth, we can choose k to be any positive integer: in particular we choose k = 
kg where 



fc = 



/l / d 2 T(t + At,t )At\ 
2 l0g25 / 3 l e 



(73) 



Equation ([73]) implies that k £ (d 2 T{t + At, t ) Ai/e)° (1) and (25/3) fe ° G (d 2 T{t + At, t ) Ai/e) o(1) because the 
square-root of a logarithm grows slower than a logarithm. We also have 

(d 2 f (t + At,t )At/e) 1/2ko G (d 2 t(t + At,i )At/e)°« (74) 

because fco increases with T(to + At,to)At/e. We then obtain the scaling given in Eq. (Til?)) by substituting these 
expressions into (|48p and using Lemma [5] □ 

For Corollary [3J it was required that there exists K > such that T'(t) < K 2 T 2 (t), because this was part of the 
conditions of Theorem [5J However, it can be expected that this requirement is automatically satisfied when {H^} is 
T-oo-pointwise-smooth, provided T is chosen appropriately (see Appendix [Bjl . 

Before concluding, we should also estimate the space-complexity of our algorithm. It follows from an analysis of 
the BACS decomposition algorithm, that the number of qubits needed for our simulation is 0(n(log* n) 2 ). Unlike the 
BACS simulation algorithm [8], this number does not depend on and e. 
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IX. CONCLUSIONS 



We introduce in this paper a pair of quantum algorithms that can be used to simulate time-dependent Hamiltonian 
evolution on quantum computers, using constant-size or adaptively chosen time steps. The adaptive time step method 
can provide superior performance, but requires that upper bounds on the norm of the Hamiltonian and its derivatives 
are known throughout the simulation. In both cases, these simulation algorithms can be performed with similar query 
complexity to the BACS simulation algorithm Q if H(t) is a sum of sufficiently smooth terms. 

We also show how to resolve pathological examples, such as our earlier example [12j wherein the higher-order 
derivatives of H(t) diverge at one point, although this method cannot be used to attain near-linear scaling unless H(t) 
is a sum of terms that are piecewise sufficiently smooth. Furthermore, we have shown that the number of operations 
used in a simulation of time-dependent Hamiltonian evolution can be reduced by using lower-order Lie- Trotter-Suzuki 
formula? to approximate time-evolution near singularities in the Hamiltonian, and higher-order formulae farther away 
from the singularities. This approach may also be useful in approximating the solutions to singular differential 
equations on classical computers. 

It may be difficult to compute some of these quantities, such as A, for some Hamiltonians. In such circumstances 
our simulation schemes can still be used, but the output state of the simulation may not be correct within an error 
tolerance of e. We recommend that heuristic testing be used to estimate whether the error is within this tolerance 
in such circumstances. Such a method could involve performing the swap test [24| between the output states of two 
separate simulations that employ distinct values of the uncertain parameter. However values such as n and the upper 
bound for ||i^|| that the oracle uses to encode the matrix elements must be known in order to perform the simulation. 

We quantify the computational complexity by the number of calls to the oracles QCol, QMatrixVal, and {T^}. 
These oracles would not typically be fundamental operations, but would be quantum subroutines consisting of se- 
quences of fundamental quantum operations. The oracles can be implemented efficiently by the quantum computer 
if each H^(t) is row-computable at every time during the simulation, and there are efficient quantum circuits for the 
basis transformations T^. 

This work leaves open several interesting avenues for investigation. One such avenue is to address the question of 
whether or not strictly linear-time quantum simulation algorithms are possible if the Hamiltonian is time-dependent 
and sufficiently smooth. In addition, it would be interesting to determine whether or not it is possible to devise an 
algorithm for which the query complexity scales poly-logarithmically with the reciprocal of the error tolerance, rather 
than sub-polynomially as our algorithm does. 



Appendix A: Proof of Lemma [T] 



In this section we present a proof of Lemma [TJ 
discretization error is bounded above by e/2. 



It provides values for both n# and n t that ensure that the 



Proof of Lemma[Ji Here we define a := At/2 nt , and define H(t) to be an approximation to the Hamiltonian Hit) that 
uses tih bit approximations to the matrix elements. We use H^j (t) to represent an approximation to the Hamiltonian 
j that is formed by taking tih bit approximations to each matrix element of the one-sparse H^ j. Similarly, we 
define f to be the closest mesh point to the time r that is still inside the interval I. In most cases, the nearest mesh 
point will be at most a distance of a/2 away from r, but there are cases where the distance can be up to ex. This 
occurs when r is near the boundary of one of the subintervals of X. Because the subintervals can have length no 
shorter than a, there will always be a mesh point within the subinterval, and the distance will not be greater than a. 
This implies that |r — t| < a. Then, using this notation, our goal in this proof is to show that values of n t and un- 
satisfying the inequalities in (|14p guarantee 



Wcxp Wexp 

]J Tlexp[-iH Mp (r p )At p }T^ p - I* exp[-i^ p , ip (f p )At p ]T M!J 
p— 1 p— 1 

where j denotes an element from a particular sequence of one-sparse Hamiltonians. 
Using Eq. (|38[) we find that, for Hermitian operators A and B, 

|i e -iA_ e -ifl|| = lim ||exp(-iA/n) n -exp(-iB/n) n \\ 

< lim n ||exp(— iA/n) — exp(— iB/n)\\ 

n— >oo 

= lim [\\A - B\\ + 0(l/n)] = \\A- B\\. 



<e/2, 



(Al) 



(A2) 
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Using this result with the Hermitian operators pt j p {T p ) and H i 



gives 



|| exp[-iH^ pJp (T p )At p ] - exp[-i J ff AlpJp (f p )Ai p ]|| < \\H^ p>jp (T p ) - i^ pJp (f p )||Ai p . 
Then we obtain, using Eq. ([58|. 

JV elt p JVcxp 

II T l e M-iH, pJp (T p )At p }T^ p - H rt p exp[-i^, ip (f p )At p ]T Mp 



(A3) 



JV r .- 



p— 1 p— 1 



(A4) 



Our next step is to bound the sum ^2 p =i |Aip|. To do so we note that, for the simulation, the time is broken up 
into r short intervals, and on each of these the Lie- Trotter-Suzuki formula Uk is used. We denote the r intervals by 
2i,22, • • • ,Ir, and the durations of these intervals by Ti, T2, . . . , T r . Using Eq. (A. 3) of Ref. [IH, if the p th exponential 
is part of the q th interval, then the duration of that exponential, At p , is at most 



|At„| < (2k/3 k )T q . 



(A5) 



The Lie- Trotter-Suzuki formula Uk is composed of 12Md 2 5 k 1 exponentials, each with a duration that is bounded 
above by (|A5|) . In addition, J2q=i Tq < At, so the total duration of the exponentials used to simulate the evolution 



in the interval I q is at most 



\ At p\ < 8A;Md 2 (5/3)' ! - 1 r, < 8kMd?(5/3) k - 1 At. 



(A6) 



Lemma [T] can be used generally in this work, because this relation does not require that the r intervals have the same 
duration, or that I is a continuous time interval. The only requirement we have used is that the Lie- Trotter-Suzuki 
integrator Uk has been used. This is to obtain the relation (|A5I) and the number of exponentials in the integrator. 
Next, using the triangle inequality, we have 



Using Taylor's theorem, we find that an upper bound for the error due to the time discretization is 

(f p )\\ < max || d t Hp(t) || a. 

By using a value of n* that satisfies (fl4| . we obtain 



H, 



mp ,jp 



(t p ) - H Ppdp (f p )\ 



< 



(32fcMd 2 )(5/3) fe - 1 Ai 



(A7) 
(A8) 

(A9) 



Next we consider the error due to the discretization of H. The matrix elements are encoded in polar form, so errors 
emerge because of inaccuracies in the modulus as well as the phase. Using the triangle inequality, the error is bounded 
above by e p + e0A max , where is the discretization error in the phase, e p is the error in the modulus, and A max is an 
upper bound for the magnitudes of the matrix elements. We choose the same number of bits to encode the modulus 
and the phase. Then taking tih to exceed the value in (|14p . the error in the modulus and phase satisfy 



e P < A 



/ (2"" /2 ) < I 
H < 2n/ (2 nH ' 2 



1 

< - 



8 (32fcMd 2 )(5/3) fe - 1 At 
2?re 



(32fcMd 2 )(5/3) fe - 1 A max At 



Using these relations, we obtain 



< £p + e0A max < 



(32fcMd 2 )(5/3) fe - 1 Ar 



Inserting Eqs. (|A6|) and (|A7j) into Eq. (JA4J) gives the error bound as 



max 
p 



^CpJp( T p) ^lUpJpfe) 



E i A *pi 



< 



2c 



x 8kMd?(5/3) k - 1 At = e/2 



P =i 



(32fcMd 2 )(5/3) fe - 1 Ai 



(A10) 
(All) 

(A12) 
(A13) 



□ 
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Appendix B: Bounds On Derivatives of T 

In our adaptive simulation method we have required that T is chosen such that there exists a constant K such 
that |9 t T(t)| < K 2 [T(t)] 2 for all times in the interval. We show in this appendix that this requirement is natural by 
demonstrating that it naturally emerges for Hamiltonians where T is chosen to be the smallest permissible function. 

Now define the set of functions {Tp} to be the smallest possible functions such that {H^} is Yp-P-pointwise- 
smooth. Taking the derivative of Y P (i) gives 



T'p(t) = lim 



st^o 



?p(t + 5t)-T P (t) 
St 



<sup{(d u \\H(u)\\)\ u = t ,-.-, (a.Har^HUVp)!^,... (duKH(u)\\ 1/{P+1) )U 



There are now two possible cases, either 



< sup < - 



i/p-i 



H ( P\t) 



P =L. 



U—t 

P+l 



(Bl) 



or 



In the first case we obtain 



and in the second case 



i/p 



< 



i/p 



> 



H (p) (u 



V(p+i) 



VO+i) 



»(*) 



i/p-i 



i/p-i 



H<*)(t) 



H^\t) 



< 



< - 

P 



2/(p+l) 



2/p 



(B2) 



(B3) 



From the definition of T P (t), \\H^(t)\\ 1 ^ p+ ^ < T P (t) for p = 0,1,..., P. However, because ||ij( p+1 ) (i)|| V(^+2) 
not included in the definition of Tp(t), we use ||ij( p+1 ^(t)|| 1 /( p+2 ' < Tp+i(f) to bound it. Then, using the fact 
that Tp+i(f) > Tp(t) and p > 1, the derivative of Tp(i) is bounded above by 



T'p(i) < [T P+1 (t)] 2 . 

A lower bound for the derivative can be obtained in the same way, giving the general result 

|T' P (t)| < [Tp +1 (t)] 2 . 



(B4) 



(B5) 



If there exists a constant K such that for all t € [to, to + At], Yp+i(f) < KTp(t), then we obtain the restriction in 
Theoreml |T' P (t)| < K 2 [T P (t)} 2 . 

In the case where T m is taken to be the smallest possible function such that {H^} is Too-oo-pointwise-smooth, this 
restriction need not be made. We see from taking the limit as P — > oo of (|B5[) that 



\?Ut)\ < { T °o(t)f 



(B6) 



This means that, if {H^} is T-oo-pointwise-smooth, the condition |T'(i)| < [T(t)} 2 should hold if T(t) is chosen 
appropriately. (It does not imply this condition, because T(t) could be chosen poorly.) 
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